Kinetics of cell division in epidermal maintenance 
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The rules governing cell division and differentiation are central to understanding the mechanisms 
of development, aging and cancer. By utilising inducible genetic labelling, recent studies have 
shown that the clonal population in transgenic mouse epidermis can be tracked in vivo. Drawing 
on these results, we explain how clonal fate data may be used to infer the rules of cell division and 
differentiation underlying the maintenance of adult murine tail-skin. We show that the rates of cell 
division and differentiation may be evaluated by considering the long-time and short-time clone fate 
data, and that the data is consistent with cells dividing independently rather than synchronously. 
Motivated by these findings, we consider a mechanism for cancer onset based closely on the model 
for normal adult skin. By analysing the expected changes to clonal fate in cancer emerging from 
a simple two-stage mutation, we propose that clonal fate data may provide a novel method for 
studying the earliest stages of the disease. 

PACS numbers: 87.17.Ee, 87.23. Cc 



I. INTRODUCTION 

A major challenge in biology is to determine how pro- 
liferating cells in developing and adult tissues behave in 
vivo. A powerful technique in solving this problem is 
clonal analysis, the labelling of a sample of cells within 
the tissue to enable their fate and that of their progeny to 
be tracked [l| . This approach gives access to information 
on proliferation, migration, differentiation (into other cell 
types) , and cell death (apoptosis) of the labelled cell pop- 
ulation. The most reliable method of labelling is through 
genetic modification leading to the expression of a re- 
porter gene in a random sample of cells. Recently it has 
become possible to activate genetic labelling at a defined 
time in transgenic mice, enabling the kinetics of labelled 
cells to be studied with single-cell resolution in vivo ||. 
From a theoretical perspective, the analysis of clonal fate 
data presents a challenging "inverse problem" in popula- 
tion dynamics: While it is straightforward to predict the 
time-evolution of a population distribution according to 
a set of growth rules, the analysis of the inverse problem 
is more challenging, open to ambiguity and potential mis- 
interpretation. 

These principles are exemplified by the mechanism of 
murine epidermal homeostasis: Mammalian epidermis is 
organised into hair follicles interspersed with interfollicu- 
lar epidermis (IFE), which consists of layers of specialised 
cells known as keratinocytes Q (see Fig. Ufa)). Prolif- 
erating cells are confined to the basal epidermal layer. 
As they differentiate into specialised skin cells, the basal 
cells withdraw from the cycle of cell proliferation and 
then leave the basal layer, migrating towards the epi- 
dermal surface from which they are ultimately shed. To 
maintain the integrity of the tissue, new cells must be 
generated to replace those lost through shedding. For 
many years, it has been thought that interfollicular epi- 
dermis is maintained by two distinct progenitor cell pop- 
ulations in the basal layer. These comprise long-lived 
stem cells (S) with the capacity to self-renew, and their 



progeny, known as transit-amplifying cells (TA), which 
go on to differentiate and exit the basal layer after sev- 
eral rounds of cell division [4(. Stem cells are also found 
in the hair follicles, but whilst they have the potential to 
generate epidermis in circumstances such as wounding, 
they do not appear to contribute to maintaining normal 
epidermis 0, @. 

The prevailing model of interfollicular homeostasis 
posits that the tissue is organised into regularly sized 
"epidermal proliferative units" or EPUs, in which a cen- 
tral stem cell supports a surrounding, clonal, popula- 
tion of transit amplifying cells, which in turn generate 
a column of overlying differentiated cells @, [1] • Several 
experimental approaches have been used to attempt to 
demonstrate the existence of EPUs, but conclusive evi- 
dence for their existence is lacking. The EPU model pre- 
dicts that slowly-cycling stem cells should be found in 
a patterned array in the IFE; cell labelling studies have 
failed to demonstrate such a pattern [9J. In chimaeric 
mice the EPU model predicts that the boundaries of mo- 
saicism in the IFE should run along the boundaries of 
EPUs; instead boundaries were found to be highly irreg- 
ular [lCj . Genetic labelling studies using viral infection 
or mutation to activate expression of a reporter gene in 
epidermal cells have demonstrated the existence of long- 
lived, cohesive clusters of labelled cells in the epidermis, 
but these clusters do not conform to the predicted size 
distribution of the EPU @, i, [Q, GJ, [H, Gj . 

Thus, until recently the means by which homeostasis 
of IFE was achieved has been unclear. However, by ex- 
ploiting inducible genetic labelling, recent studies have 
allowed the fate of a representative sample of progenitor 
cells and their progeny to tracked in vivo hj. As well 
as undermining the basis of the stem/TA cell hypothesis, 
the range of clone fate data provide the means to infer 
the true mechanism of epidermal homeostasis. In partic- 
ular, these investigations indicate that the maintenance 
of IFE in the adult system conforms to a remarkably 
simple birth-death process involving a single progenitor 
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FIG. 1: (Color online) (a) Schematic cross-section of murine 
interfollicular epidermis (IFE) showing the organisation of 
cells within different layers and indicating the architecture 
of typical labelled clones. Proliferating cells (grey) are con- 
fined to the basal layer (labelled i); differentiated cells mi- 
grate through the superbasal layers (ii), where they flatten 
into cornified cells, losing their nuclei and assembling a corni- 
fied envelope (green) (iii), eventually becoming shed at the 
surface. The shaded regions (yellow) indicate two distinct 
clones, the progeny of single basal layer cells labelled at in- 
duction. While the clone on the right retains at least one 
labelled cell in the basal layer, the clone on the left hand side 
has detached from the basal layer indicating that all of the 
cells have stopped proliferating. The former are designated as 
"persisting clones" and contribute to the clone size distribu- 
tions, while the latter, being difficult to resolve reliably, are 
excluded from experimental consideration, (b) Typical exam- 
ple of a clone acquired at a late time point, viewed from the 
basal layer surface. Cell nuclei are labelled blue; the heredi- 
tary clone marker (EYFP) appears yellow. Scale bar: 20/^m 



cell compartment. Expanding upon the preliminary the- 
oretical findings of Ref . [HI , the aim of this paper is to 
elucidate in full the evidence for, and the properties of, 
the model of epidermal maintenance, and to describe the 
potential of the system as a method to explore early sig- 
natures of carcinogenic mutations. 



A. Background: Experimental Methodology 

To organise our discussion, we begin with an overview 
of the experimental arrangement, referring to Ref. fl5j 
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FIG. 2: (Color online) Top: Theoretical lineage for the first 
12 weeks post-labelling of (a) a detached clone in which all 
cells have undergone a transition to terminal differentiation 
by week 12, and (b) a persisting clone in which some of the 
cells maintain a proliferative capacity, according to model ([2} . 
Circles indicate progenitor cells (P), differentiated cells (D), 
and suprabasal cells (SB). Note that, because the birth-death 
process ([2]) is Markovian, the lifetime of cells is drawn from 
a Poisson distribution with no strict minimum or maximum 
lifetime. The statistics of such lineage trees do not change 
significantly when we account for a latency period between 
divisions that is much shorter than the mean cell lifetime (see 
discussion in section III C[) . Bottom: The total number of 
proliferating, differentiated and supra-basal cells for the two 
clones as a function of time. 



for technical details of the experimental system. To gen- 
erate data on the fate of individual labelled cells and 
their progeny, hereafter referred to as clonal fate data, 
inducible genetic marking was used to label a sample of 
cells and their progeny in the epidermis of transgenic 
mice. The enhanced Yellow Fluorescent Protein (EYFP) 
label was then detected by confocal microscopy, which 
enables 3D imaging of entire sheets of epidermis. Low- 
frequency labelling of approximately 1 in 600 basal-layer 
epidermal cells at a defined time was achieved by using 
two drugs to mediate a genetic event which resulted in 
expression of the EYFP gene in a cohort of mice. This 
low efficiency labelling ensures that clones are unlikely 
to merge (see discussion in section III A[) . By analysing 
samples of mice at different time points it was possible 
to analyse the fate of labelled clones at single cell reso- 
lution in vivo for times up to one year post-labelling in 
the epidermis (see, for example, Fig. [JJb)) 

With the gradual accummulation of EYFP levels, the 
early time data (less than two weeks) reveals a small in- 
crease in the number of labelled clones containing one 
or two cells. At longer times, clones increase in size 
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FIG. 3: (Color online) Mean number of basal layer cells in 
persisting clones. The experimental data (circles) show an 
inexorable increase in the size of persisting clones over the 
entire time course of the experiment. The behaviour at short 
times (from 2 — 6 weeks) and at long times (beyond 13 weeks) 
follows the two simple analytical approximations described in 
the main text (lower and upper dashed curves). For times ear- 
lier than two weeks (referring to section HlD[l . clones remain 
approximately one cell in size. The experimental data are 
consistent with the behaviour predicted by process @ (black 
line) when it is assumed that only A-type cells are labelled 
at induction. In contrast, assuming that A and B type cells 
label in proportion to their steady-state population leads to 
an underestimate of average clone size between two and six 
weeks (lower curve, red online), as does the assumption that 
type B cells label with better efficiency (not shown). Inset: 
The underlying distribution of basal cells per clone at 2 weeks 
and 26 weeks post-labelling. The data is binned by cell count 
in increasing powers of 2. 



sisting clone population. This result is incompatible with 
any model in which the IFE is supported by a pop- 
ulation of long-lived stem cells. With the latter, one 
would expect the number density of persisting clones to 
reach a non-zero minimum (commensurate with the la- 
belling frequency of stem cells) while the average clone 
size would asymptote to a constant value characteristic 
of a single epidermal proliferative unit. We are therefore 
lead to abandon, or at least substantially revise, the or- 
thodox stem/TA cell hypothesis and look for a different 
paradigm for epidermal maintenance. 

But, to what extent are the clone fate data amenable 
to theoretical analysis? Indeed, the application of pop- 
ulation dynamics to the problem of cell kinetics has a 
long history (see, e.g., with studies 

of epidermal cell proliferation addressed in several pa- 
pers [2fl l2ll . I22I |23| . However, even in the adult sys- 
tem, where cell kinetics may be expected to conform to a 
"steady-state" behaviour, it is far from clear whether the 
cell dynamics can be modelled as a simple stochastic pro- 
cess. Regulation due to environmental conditions could 
lead to a highly nonlinear or even non-local dependence 
of cell division rates. Indeed, a priori, it is far from clear 
whether the cell kinetics can be considered as Markovian, 
i.e. that cell division is both random and independent of 
the past history of the cell. Therefore, instead of trying 
to formulate a complex theory of cell division, taking ac- 
count of the potential underlying biochemical pathways 
and regulation networks [20| . we will follow a different 
strategy looking for signatures of steady-state behaviour 
in the experimental data and evidence for a simple under- 
lying mechanism for cell fate. Intriguingly, such evidence 
is to be found in the scaling properties of the clone size 
distribution [Tj|. 



while cells within clones begin to migrate through the 
suprabasal layers forming relatively cohesive irregular 
columns (see Fig.QJa)). 

The loss of nuclei in the cornified layer (fig. [T]) makes 
determination of the number of cornified layer cells in 
larger clones by microscopy unreliable. Therefore, to 
identify a manageable population, attention was focused 
on the population of basal cells in "persisting clones", 
defined as those labelled clones which retain at least one 
basal layer cell, such as is exemplified in the theoreti- 
cal lineage maps in Fig. [2j After two weeks, the density 
of persisting clones was seen to decrease monotonically 
indicating that the entire cell population within such 
clones had become differentiated and the clone detached 
from the basal layer (shown schematically in figs. []Ja) 
and[2{a)). However, the population of persisting clones 
showed a steady increase in size throughout the entire 
duration of the experiment. 

To what extent are the clone fate data consistent with 
the orthodox stem/TA cell model of epidermal mainte- 
nance? Referring to Fig. [3J one observes an inexorable 
increase in the average size of an ever-diminishing per- 



B. Scaling 

To identify scaling characteristics, it is necessary to 
focus on the basal layer clone size distribution, P n (t), 
which describes the probability that a labelled progen- 
itor cell develops into a clone with a total of n basal 
layer cells at a time t after labelling. (Note that, in gen- 
eral, the total number of cells in the supra-basal layers 
of a clone may greatly exceed the number of basal layer 
cells.) With this definition, Po(t) describes the "extinc- 
tion" probability of a clone, i.e. the probability that all 
of the cells within a labelled clone have migrated into the 
supra-basal layers. To make contact with the experimen- 
tal data, it is necessary to eliminate from the statistical 
ensemble the extinct clone population (which are difficult 
to monitor experimentally) and single-cell clones (whose 
contribution to the total ensemble is compromised by the 
seemingly unknown relative labelling efficiency of prolif- 
erating and post-mitotic cells at induction), leading to a 
reduced distribution for "persisting" clones, 

ppcrs./.x _ Pn(t) 

n*> 2 w- 1 _ Po{t) _p lit y 
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Then, to consolidate the data and minimise fluctuations 
due to counting statistics, it is further convenient to bin 
the distribution in increasing powers of 2, 

2 k 

n=2 k ~ 1 +l 

i.e. Vf eTS '(t) describes the probability of having two cells 
per clone, "P^ ers ' (t) describes the probability of having 
3-4 cells per clone, and so on. Referring to Fig. [U one 
may see that, after an initial transient behaviour, the 
clone size distribution asymptotes in time to the simple 
scaling form, 

rr s -(t) = f(2 k /t). (i) 

This striking observation brings with it a number of 
important consequences: As well as reinforcing the inap- 
plicability of the stem cell/TA cell hypothesis, such be- 
haviour suggests that epidermal maintenance must con- 
form to a simple model of cell division. The absence 
of further characteristic time-scales, beyond that of an 
overall proliferation rate, motivates the consideration of 
a simple kinetics in which only one process dictates the 
long-time characteristics of clonal evolution. 

Moreover, from the scaling observation one can also 
deduce two additional constraints: Firstly, in the long- 
time limit, the average number of basal layer cells within 
a persisting clone increases linearly with time, viz. 

(n) pcrs . = E" P ™>TW- / dn-f(n/t) ext. 

n>2 J ° 1 

Secondly, if we assume that labelled progenitor cells are 
representative of all progenitor cells in the epidermis, 
and that the population of clones with only one basal 
layer cell is not "extensive" (i.e. Hindoo P\(t) — 0), this 
means that, in the long-time limit, the clone persistence 
probability must scale as 1 — Pq (t) oc 1/t such that 

n 

where the constant, p, is given by the fraction of prolif- 
erating cells in the basal layer. Without this condition, 
one is lead to conclude that the labelled population of 
basal layer cells either grows or diminishes, a behaviour 
incompatible with the (observed) steady-state character 
of the adult system. 

Although the manifestation of scaling behaviour in the 
clone size distributions gives some confidence that the 
mechanism of cell fate in IFE conforms to a simple non- 
equilibrium process, it is nevertheless possible to conceive 
of complicated, multi-component, models which could 
asymptote to the same long-time evolution. To further 
constrain the possible theories, it is helpful to draw on 




FIG. 4: (Color online) Time dependence of the grouped size 
distribution of persisting clones, V^ CIS '(t), plotted as a func- 
tion of the rescaled time coordinate t/2 k i— » t. The data points 
show measurements (extracted from data such as shown in 
Fig- EJinset), given fully in ref. [H|]), while the solid curves 
show the probability distributions associated with the non- 
equilibrium process for the basal-layer clone population as 
obtained by a numerical solution of the Master equation (3). 
(Error bars refer to standard error of the mean). At long 
times, the data converge onto a universal curve (dashed line), 
which one may identify with the form given in eq. [7] The 
rescaling compresses the time axis for larger clones, so that 
the large-clone distributions appear to converge much earlier 
onto the universal curve. 



additional experimental observations (l5| : Firstly, im- 
munostaining of clones with a total of two cells (using 
the proliferation marker Ki67 and, separately, the repli- 
cation licensing factor cdc6) reveals that a single cell 
division may generate cither one proliferating and one 
non-proliferating daughter through asymmetric division, 
or two proliferating daughters, or two non-proliferating 
daughters (cf. [24j, |2fj|, |26[ ) . Secondly, three-dimensional 
imaging of the epidermis reveals that only 3% of mitotic 
spindles lie perpendicular to the basal layer indicating 
that divisions may be considered to be confined to the 
basal layer, confirming the results of earlier work that 
indicates a dividing basal cell generates two basal layer 
cells [Hj]. 

This completes our preliminary discussion of the exper- 
imental background and phenomenology. In summary, 
the clone fate data reveal a behaviour wholely incompat- 
ible with any model based on the concept of long-lived 
self-renewing stem cells. The observation of long-time 
scaling behaviour motivates the consideration of a simple 
model based on a stochastic non-equilibrium process and 
is indicative of the labelled cells being both a representa- 
tive (i.e. self-sustaining) population and in steady-state. 
In the following, we will develop a theory of epidermal 
maintenance which encompasses all of these observations. 
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II. THEORY OF EPIDERMAL MAINTENANCE 
A. Model 

Taken together, the range of clonal fate data and the 
observation of symmetric and asymmetric division are 
consistent with a remarkably simple model of epidermal 
homeostasis involving only one proliferating cell compart- 
ment and engaging just three adjustable parameters: the 
overall cell division rate, A; the proportion of cell divi- 
sions that are symmetric, r; and the rate of transfer, T, of 
non-proliferating cells from the basal to the supra-basal 
layers. To maintain the total proliferating cell popula- 
tion, a constraint imposed by the steady-state assump- 
tion, we have used the fact that the division rates asso- 
ciated with the two channels of symmetric cell division 
must be equal. Denoting the proliferating cells as type 
A, differentiated basal layer cells as type B, and supra- 
basal layer cells as type C, the model describes the non- 
equilibrium process, 




A Prob. r 

B Prob. 1 - 2r 

B Prob. r 



(2) 



Finally, the experimental observation that the total basal 
layer cell density remains approximately constant over 
the time course of the experiment leads to the additional 
constraint that 



r = 



1-P 



A. 



reducing the number of adjustable parameters to just 
two. 

By ignoring processes involving the shedding of cells 
from the surface of the epidermis, the applicability of the 
model to the consideration of the total clone size distri- 
bution is limited to appropriately short time scales (up 
to six weeks post-labelling). However, if we focus only 
on the clone size distribution associated with those cells 
which occupy the basal layer, the model can be applied 
up to arbitrary times. In this case, the transfer process 

must be replaced by one in which B — > 0. In cither case, 
if we treat all instances of cell division and cell transfer 
as independent stochastic events, a point that we shall 
revisit later, then the time evolution associated with the 
process ^ can be cast in the form of a Master equation. 
Defining P nA ,n B {t) as the probability of finding n A type 
A cells and «b type B cells in a given clone after some 
time t, the probability distribution evolves according to 
the Master equation: 



dtPn A ,n B ~ r ^ i( n A ~ l)-fnA — 1|71B — n APn A ,n B ] 

+r\[(n A + l)P„ A+l! „ B _ 2 - n A P nA , nB 
+ (1 - 2r)X[n A P nA 

,71b — 1 - n A P nA 

,«bJ 

+r[(n B + l)P„ A ,„ B+ i - n B P nA ,n B ] ■ 



(3) 



If we suppose that the basal layer cells label in proportion 
to their population, the latter must be solved subject to 
the boundary condition P nA ,n B {0) — pSn A ,i$n B ,a + (1 - 
p)^n A ,o^n B ,i- Later, in section [II Dl we will argue that 
the clone size distribution is compatible with a labelling 
efficiency which favours A over B type cells. Either way, 
by excluding single cell clones from the distribution, this 
source of ambiguity may be safely eliminated. Although 
the Master equation (and its total cell number general- 
isation) is not amenable to exact analytic solution, its 
properties can be inferred from the consideration of the 
A cell population alone for which an explicit solution may 
be derived. 

When considered alone, A type cells conform to a sim- 
ple set of rate laws, 



2rX 



A + A 



Prob. 
Prob. 



1/2. 
1/2. 



(4) 



an example of a Galton- Watson process, long known 
to statisticians (see, e.g., Ref. |27|). In this case, the 
probability distribution, which is related to that of the 
two-component model through the relation p nA (t) = 
Y^L=a Pn A ,n B ('t)> can be solved analytically. (Here, we 
have used a lower case p to discriminate the probability 
distribution from its two-component counterpart.) For 
an initial distribution p„ A (0) = 5n A ,i it m &y be shown 
that [H, 



Pn A (t) = 1 



1 

rXt 



-(riA+l) 



1 



UA = , 

7 n A > . 



(5) 



From this system and its associated dynamics, one can 
draw several key implications: 



1. Epidermis is maintained through an ever-decreasing 
clonal population: 

Starting with a single labelled cell, the Galton- Watson 
process predicts that the persistance probability of the 
resulting clone (i.e., in this case, the probability that the 
clone retains at least one proliferating cell), is given by 



Pn A >o = 1 ~Po(t) = 



1 



1+rXt 



i.e. as with the experiment, the persistance probability 
of a clone decays monotonically, asymptoting to the form 
1 — po(t) oc 1/t at time scales t 3> 1/rA, the time scale for 
symmetric division. Applied to the experimental system, 
this suggests that labelled clones continue to detach from 
the basal layer indefinitely. At the same time, defining 



pcrs. /,\ 

<a>o(*) = 



Pn A (t) 
1 -«>(*) 



as the size distribution of persisting clones, the mean 
number of basal layer cells in a persisting clone grows 
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steadily as 



and so 



("a) 



pcrs. 



Xt. 



such that the overall cell population remains constant, 
viz. (tia) = J2 < ^=o riA P n A(t) — I; i- e - the continual ex- 
tinction of clones is compensated by the steady growth 
of persisting clones such that the average number of pro- 
liferating cells remains constant: given enough time, all 
cells would derive from the same common ancestor, the 
hallmark of the Galton- Watson process [38[ . 

This linear increase in clone size may lead one to worry 
about neighbouring clones coalescing. Fortunately, the 
continual extinction of clones ensures that the fraction 
of clones conjoined with their neighbours remains small 
and of same order as the initial labelling density [39| . 
The fact that this fraction is constant is again indicative 
of the steady-state condition maintained throughout the 
experiment. 



2. Larger clones begin to exhibit the stability of the 
macroscopic system: 

If, at some instant, a clone is seen to have, say, Na 
proliferating cells then, after a further time t, its size will 
fluctuate as 



((n A -(»A)) 2 ) 1/2 
(«a) 



2r\t 

~/v7 



Thus clones (as defined by the A cell population) will 
maintain an approximately stable number of cells provid- 
ing t <c Na/t\. For larger clones this time may exceed 
the lifetime of the system. At the limit where macro- 
scopic sections of the basal layer are considered, the sta- 
tistical fluctuations are small. The increased stability 
of larger clones also explains the surprising prediction 
that, given enough time, all clones eventually become 
extinct (viz. limt_ >0 oPn>o(^) — 0). Calculated explic- 
ity, the extinction probability for a clone of size Na ~S> 1 
scales as po(t) ~ e - N A/r\t approaching unity at long 
times. However, because this extinction probability is 
small when t <C Na/tX, a large enough clone may easily 
persist beyond the lifetime of the system. 



lim 

t»l/rA 



exp 



2r\t 



-exp 



rXt 



(7) 



i.e. the probability distribution acquires the scaling form 
found empirically. Referring to Eq. (TI]) , we can therefore 
deduce the form of the scaling function, 



f(x) = exp[— px/2rX] — exp[— px/rX] 



(8) 



As a result, at long times, the average basal layer pop- 
ulation of persisting clones becomes proportional to the 
average number of proliferating cells per clone, (n) pors . = 
(1 + r\t)/p, a behaviour consistent with that seen in ex- 
periment (see Fig. [3]). 



4- The creation and transfer of differentiated cells dictates 
the short-time behaviour of the clone size distribution: 

In fitting the model to the data (see below), we will 
find that the rates A and T at which differentiated cells 
are created and then transferred into the super-basal re- 
gion are significantly larger than the rate of symmetric 
division rX, which dictates the long-time behaviour of 
the clone size distribution. In this case, at early times 
(t ^ the clone size distributions are dominated 

by the differentiation and transfer rates, which remain 
prominent until the population of labelled differentiated 
cells associated with each proliferating cell reaches its 
steady-state value of (1 — p)/p. One may therefore infer 
that, at short times, the mean number of basal layer cells 
in clones arising from proliferating cells is given by 



lim (n)r 
t«i/r 



l/p-(l/p-l)e 



n 



and that the early-time clone size distribution is Poisson- 
distributed, viz. 



lim P£™"(t) 
t<i/r n - 



((»> 



pcrs. 



1) (n-1)! 



(9) 



B. Fit to the data 



3. The properties of the proliferating cell population 
dictates the behaviour of the entire clone size distribution: 

At asymptotically long times, one may show [4(| that 
the full probability distribution for finding n = ua + 
ris cells within a persisting clone scales in proportion to 
pP- s -(t), viz. 



lim . 

t»l/rA 



rXt 



exp 



pn 
rXt 



(6) 



With these insights it is now possible to attempt a fit of 
the model to the data. Referring to Fig. one may infer 
the rate of cell division A from the short-time data, and 
the symmetric division rate r X from the long-time scaling 
data. In particular, taking the fraction of proliferating 
cells in the basal layer to be p = 0.22, a figure obtained 
experimentally by immunostaining using Ki67 }15j . a fit 
of Eq. ([9]) to the short-time data (fig.OJa)) is consistent 
with a transfer rate of T = 0.31/week which, in turn, im- 
plies a rate of cell division of A = 1.1/week. Furthermore, 
by plotting the long-time, large-fc, size-distributions in 
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FIG. 5: (Color online) (a) Fit of Eq. Q to the short-time clone size distributions. At early times (a), the data is optimally 
fitted by Eq. © using the value A = f.f/week and the empirical value p = 0.22 (solid lines show fit). To ensure integrity 
of the analysis, data for times earlier than week 2 have been excluded (see section III D[l . (b) Linearisation of the long-time 
asymptotic data using the "inverse" scaling function f~ 1 (2 k /t) for t > 13 weeks and k > 3 (see main text), (c) Likelihood of 
the overall division rate A (red) and the symmetric division rate r\ (black), as assessed from a x 2 test of the numerical solution 
to Eq. © 29]. A fit to the basal-layer clone size distribution alone (dashed) is less discriminatory than a simultaneous fit to 
both the basal- layer and total clone size distributions (solid curves). The likelihood of rA is shown for the optimal value of 
A, and vice- versa. Inset: Referring to section Til CI the likelihood is plotted against the duration of a latency period (r m i n .) 
immediately following cell division, and assuming that division events are otherwise independent (see main text). 



terms of the "inverse" to the scaling function, 
f-\2 k /t) = (21n[(l-(l-/(2 fc A)) 1/2 )/2]) _1 
"(l-(l-Pr S -W) 1/2 )/2 
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the data converge onto a linear plot (Fig. Efb)). The 
resulting slope takes the value —rX/p, from which we may 
infer the symmetric division rate rA = 0.09 ± 0.01/week, 
and r = 0.08 ±0.01. 

These figures compare well with an optimal fit of 
the entire basal layer clone size distribution (Fig. Q|, 
obtained by numerically integrating the Master equa- 
tion ([3]). The fitting procedure is shown in Fig. [5^c) 
(solid curves), where the likelihood of the model is eval- 
uated for a range of values of A and rA, as assessed from 
a x 2 test of the model solution 29]. One may see that 
the likelihood is maximised with an overall division rate 
of A = 1.1 /week and a symmetric division rate in the 
range rA = 0.1 ± 0.01/week, thus confirming the validity 
of the asymptotic fits. Moreover, the corresponding fit of 
both the basal layer distribution and the total clone size 
distribution, including both basal and supra-basal cells, 
is equally favourable (Fig. G3c), dashed). Thus, in the 
following sections we shall use the asymptotically fitted 
value of r = 0.08, however any choice of the parameter 
in the range r = 0.08 — 0.10 gives similar results. 

Although the comparison of the experimental data 
with the model leaves little doubt in its validity, it is 



important to question how discerning is the fit. By it- 
self, the observed increase in the size of persisting clones 
is sufficient to rule out any model based on long-lived 
self-renewing stem cells, the basis of the orthodox EPU 
model. However, could one construct a more complicated 
model, which would still yield a similar fit? Certainly, 
providing the long-time evolution is controlled by a sin- 
gle rate-determining process, the incorporation of further 
short-lived proliferating cell compartments (viz. transit- 
amplifying cells) would not affect the observed long-time 
scaling behaviour. However, it seems unlikely that such 
generalisations would provide an equally good fit to the 
short-time data. 

More importantly, it is crucial to emphasize that the 
current experimental arrangement would be insensitive 
to the presence of a small, quiescent, long-lived stem cell 
population. Yet, such a population could play a crucial 
role in non-steady state dynamics such as that associated 
with wound healing or development. We are therefore 
led to conclude that the range of clone fate data for nor- 
mal adult IFE are consistent with a simple (indeed, the 
simplest) non-equilibrium process involving just a single 
progenitor cell compartment. 



C. Stochastic behaviour of cell division 

At this stage, it is useful to reflect upon the sensitivity 
of the model to the stochasticity assumption applied to 
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Cell lifetime Time (weeks) 

FIG. 6: (Color online), (a) Examples of progenitor cell cycle-time distributions with the same average cycle time 1/A, (A = 
1.1/week), and with a latency period of T m m. = 12 hours introduced between consecutive cell divisions (hashed region). The 
case k — 1 corresponds to a model of independent cell division such as assumed in section III Al but now accounting for an 
initial latency period. The case k — * oo (black dashed) corresponds to all cells having an exact cell cycle-time of 1/A. Note that 
small values of k allow for both very short and very long cycle times, (b) Using Monte-Carlo simulations of process ([2j|, the 
clone size distributions predicted by each of the different cycle-time distributions in (a) are compared with the empirical data. 
Data points show the size distribution of persisting clones including supra-basal layer cells over the first 6 weeks post-labelling 
(extracted from data given fully in Ref. [15| ; for legend see Fig. [4} , and the theoretical curves correspond to the same legend 
as in (a). All of the models give an optimal fit with the same value of A = 1.1/week, r — 0.08. 



the process of cell division. Clearly, the scaling behaviour 
(Eq. depends critically on the statistical independence 
of successive cell divisions; each cell division results in 
symmetric/asymmetric cell fate with relative probabili- 
ties as detailed in |4}. But, to what extent would the 
findings above be compromised if the cell cycle-time, i.e. 
the time between consecutive cell divisions, were not de- 
termined by an independent stochastic process? This 
question may have important ramifications, because the 
assumption of independent cell division, used in formu- 
lating the Master equation ([3]), introduces a manifestly 
unphysical behaviour by allowing cells to have arbitrar- 
ily short cycle times. Moreover, although a wide dis- 
tribution of cell cycle-times has been observed for hu- 
man keratinocytes in vitro [30], it is possible that that 
keratinocytes in vivo may divide in synchrony, giving a 
cell cycle-time distribution narrowly centered about the 
mean (1/A). In the following, we shall address both of 
these points: Firstly, we shall show that, up to some 
potential latency period (the time delay before a newly- 
divided cell is able to divide again), consecutive cell di- 
visions occur independently as an asychronous, Poisson 
process. Secondly, while the data is insufficient to detect 
a latency period of 12 hours or less between consecutive 
cell divisions, the data does discriminates against a pe- 
riod lasting longer than 24 hours. 

To investigate the degree to which the model is sen- 
sitive to the particular cell cycle-time distribution, let 
us revisit the original model of independent cell divi- 
sion with several variations: Firstly, we introduce a la- 
tency period of r m i n immediately following cell division, 



in which daughter cells cannot divide. This biologically- 
motivated constraint renders a more complicated yet 
more realistic model of cell division than the idealised 
system studied in the previous section. Motivated by 
observations of the minimal cycle-time of (human) ker- 
atinocytes [30|, where a latency period of T m ; n ~ 10 hours 
was observed in vitro, we shall here consider the a range 
latency periods of up to 48 hours. Secondly, we com- 
pare the empirical clone size distributions with a model 
where all progenitor cells have a cycle-time of exactly 
1/A, i.e. where cells within each clone divide in perfect 
synchrony. Finally, we shall investigate a range of inter- 
mediate models with different distributions of progenitor 
cell cycle-time (see Fig. [B£a)). 

Technically, the resulting clone size distributions may 
be evaluated through Monte Carlo simulations of the 
non-equilibrium process @ with the cycle-time r of each 
proliferating cell selected at random from a Gamma dis- 
tribution of the form 



where f = 1/A — r m i n . is the average time to division 
following the initial latency period r m i n ., and k is the 
"shape parameter" of the Gamma distribution. In par- 
ticular, the choice of shape parameter n = 1 corresponds 
to the exponential distribution which characterises the 
independent cell cycle-time distribution, whereas k — > oo 
describes the case in which all A-cells have an exact cycle- 
time of 1/A (see Fig.JSJa)). Then, to reflect the assump- 
tion that initially-labelled, spatially separated, progen- 
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itor cells have uncorrelated cell cycles, the time to the 
initial division event post-labelling is adjusted by a ran- 
dom time r £ [0, 1/A]. Finally, for an unbiased compar- 
ison of the models, we optimise the value of A for each 
model separately against the empirical data, whilst keep- 
ing rX = const, to ensure an optimal fit of the long-time 
data, as discussed below. 

The resulting clone size distributions are shown in 
Fig. Hlb) , where the case of independent division follow- 
ing a 12-hour latency (re = 1) and the exact cycle-time 
case (re — > oo) are compared to the empirical total clone 
size distribution, which includes both basal and supra- 
basal (type C) cells, over the first 6 weeks post-labelling. 
Two intermediate cases are also shown for comparison 
(re = 2,10). Focusing first on the results for the case 
re = 1 , which bears closest resemblance to the Markovian 
model analysed using the Master Equation Q, one may 
see by inspection that the quality of the fit to the data 
remains good even when the effects of a latency period 
between cell divisions is taken into account. More rig- 
orously, a likelihood analysis reveals that the two cases 
are statistically indistinguishable (see Fig. (He), inset), 
which indicates that the duration of a latency period of 
Tmin. < 12 hours is beyond the current empirical reso- 
lution. However, referring to Fig. (inset), a similar 
analysis of longer latency periods reveals that for periods 
of T m j n . > 24 hours, the fit to the data is significantly 
poorer. 

Turning next to the predicted basal-layer clone size 
distributions at late times (t > p/fX) (not shown), one 
may see that all of the proposed distributions asymptot- 
ically converge: Starting with exactly one cell, then the 
moment-generating function G(q, s) — Y^=oPni. s )l n as ~ 
sociated with the A cell population distribution p n (s) 
after s cell cycles satisfies the recursion relation [31j : 

G(q,S + l)-G(q,s) = r(G(q,s)-l) 2 , 

which asymptotes to the continuous master equation 
lim s2> i d s G(q, s) = r (G(q, s) — l) 2 , with the relative 
magnitude of the leading-order correction dropping off as 
1/s. But with s = At, this equation is simply the master 
equation for the moment-generating function associated 
with the original model, Eq. [5l and so the two models 
converge. One may therefore conclude that, beyond the 
first several weeks of the experiment (t ^> 1/A), the fit to 
the data is sensitive only to the average cycle time of pro- 
genitor cells. With this in mind, we note that for the case 
of perfectly synchronous cell division, an optimal (albeit 
poor) numerical fit was obtained when A = 1.2/ week, a 
figure that compares well with the fit for the indepen- 
dent case. It appears therefore that the predicted average 
cell division rate (X) is insensitive to the shape of the cell 
cycle distribution. 

Finally, let us turn to the early time behaviour (t ~ 
1/A), where the predicted distributions are distinct. Re- 
ferring to Fig. HIb), one may see, at 2-4 weeks post- 
labelling, that relatively large clones (5 — 8 cells) appear 
earlier than expected by a model assuming synchronous 



division, and that, compared with the same model, a 
sizeable proportion of small clones (e.g., 2 cells) lingers 
on for far longer than expected. The same behaviour is 
observed for the basal layer clone size distribution (not 
shown). One may therefore infer that cell division con- 
forms to a model of independent rather than synchronous 
division, allowing for some progenitor cells to divide un- 
usually early, and for others to remain quiescent for an 
unusually long period of time. 

In summary, we have established that, following di- 
vision, progenitor cells do not divide for a period that 
is likely to last up to 12 hours, and not more than 24 
hours. After this latency period, the data is consistent 
with cells switching to a mode of independent, asyn- 
chronous, cell division. These results shed light on why 
the simple model of independent cell division presented 
in section Hi Al succeeds in producing such a remarkable 
fit to the data. 



D. Labelling efficiency and EYFP accumulation in 
basal cells 

Although the integrity of the fit of the model to the 
data provides some confidence in its applicability to the 
experimental system, its viability as a model of epider- 
mal homeostasis rests on the labelled clone population 
being representative of all cells in the IFE. Already, we 
have seen that the model, and by inference, the labelled 
clone population, has the capacity to self-renew. How- 
ever, the slow accumulation of EYFP after induction, 
together with the question of the relative labelling effi- 
ciency of the two basal layer cell types, leaves open the 
question of the very short-time behaviour. Accepting the 
validity of the model, we are now in a position to address 
this regime. 

In doing so, it is particularly useful to refer to the time 
evolution of clone size as measured by the average num- 
ber of basal cells in a persisting clone. As expected from 
the scaling analysis discussed in section ITBl a compari- 
son of the experimental data with that predicted by the 
proposed cell kinetic model shows a good agreement at 
long times (Fig. [3]). However, comparison of the data at 
intermediate time-scales provides significant new insight. 
In particular, if we assume equal labelling efficiency of 
progenitor and differentiated cells, i.e. that both cell 
types label in proportion to their steady-state popula- 
tion (shown as the lower (red) curve in the Fig. [3]), then 
there is a substantial departure of the predicted curve 
from the experimental data for times of between two and 
six weeks. Intriguingly, if we assume that differentiated 
cells simply don't label, then the agreement of the data 
with theory is excellent from two weeks on! We are there- 
fore lead to conclude that, at least from two weeks, all 
labelled clones derive from progenitor cells labelled at 
induction. 

With this in mind, we may now turn to the average 
clone size as inferred from the data at two days and one 
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week. Here one finds that the model appears to sub- 
stantially over-estimate the clone size. Indeed, Fig. [3] 
suggests that the average clone size is pinned near unity 
until beyond the first week post-labelling, i.e. the rela- 
tive population of single-cell clones is significantly larger 
than expected at one week, yet falls dramatically to 
the theoretical value at two weeks. Referring again to 
the slow accumulation of EYFP, can one explain the 
over-representation of single-cell clones at one week post- 
labelling? At one week, two-cell clones are observed soon 
after cell division, and thus express lower concentrations 
of EYFP compared to single-cell clones. As a result they 
may be under-represented. At later times, all labelled 
clones become visible as EYFP concentration grows, ex- 
plaining the coincidence of experiment and theory at two 
weeks. It follows, of course, that the size distributions 
at later time points would be unaffected by slow EYFP 
accumulation. However, a full explanation of this effect 
warrants further experimental investigation, and is be- 
yond the scope of this paper. 



III. MANIFESTATION OF MUTATIONS IN 
CLONAL DISTRIBUTIONS 

Having elucidated the mechanism of normal skin main- 
tenance, it is interesting to address its potential as a 
predictive tool in clonal analysis. Conceptually, the ac- 
tion of mutations, drug treatments or other environmen- 
tal changes to the tissue can effect the non-equilibrium 
dynamics in a variety of ways: Firstly, a revision of cell 
division rates or "branching ratios" (i.e. symmetric vs. 
asymmetric) of all cells may drive the system towards 
either a new non-equilibrium steady-state or towards a 
non-steady state evolution resulting in atrofication or un- 
constrained growth of the tissue. (The development of 
closed non steady-state behaviour in the form of limit 
cycles seems infeasible in the context of cellular struc- 
tures.) Secondly, the stochastic revision of cell division 
rates or branching ratios of individual cells may lead to 
cancerous growth or extinction of a sub-population of 
clones. The former may be referred to as a "global per- 
turbation" of the cell division process while the second 
can be referred to as "local" . In both cases, one may ex- 
pect clonal analysis to provide a precise diagnostic tool 
in accessing cell kinetics. To target our discussion to the 
current experimental system, in the following we will fo- 
cus on the action of a local perturbation in the form of 
a carcinogenic mutation, reserving discussion of a global 
perturbation, and its ramifications for the study of drug 
treatment, to a separate publication. 

Let us then consider the action of a local perturba- 
tion involving the activation of a cancer gene in a small 
number of epidermal cells, which leads to the eventual 
formation of tumours. In the experimental system, one 
can envisage the treatment coinciding with label induc- 
tion, for example by simultaneously activating the EYFP 
and the cancer gene. In this case, clonal fate data should 



simply reflect a modified model of cell proliferation lead- 
ing to the eventual failure of the steady-state model of 
tissue maintenance. 



A. A simple model of carcinogenesis 

To quantify the process of cancer onset, we start by 
establishing the simplest possible changes to process ([2]) 
which may be associated with tumour growth. Cancer is 
widely held to be a disease caused by genetic instability 
that is thought to arise when a progenitor cell undergoes 
a series of mutations [H, [H, [HJ ■ As a result, cells within 
the mutant clone prefer to proliferate, on average, over 
processes leading to terminal differentiation or death. In 
this investigation we shall consider a "simple" cancer re- 
sulting from two rate-limiting mutations: Referring to 
our proposed labelling experiment, the controlled induc- 
tion of a cancer-causing mutation during label induction 
defines the first mutation; a second, rate-limiting step 
then occurs with the stochastic occurrence of a second 
cancer causing mutation. Examples of the first type of 
mutation may be genes that affect the ability of a cell to 
respond to genetic changes of the cell, e.g. p53, whilst the 
second mutation may be of a gene that affects clone fate 
such as the Ras oncogene [33[ . We may therefore distin- 
guish between "stage one" mutated cells, which maintain 
the steady-state, and "stage two" cells, which have the 
capacity for tumour formation. 

The resulting process of cell proliferation is set by three 
parameters: The overall rate of mutation v from a stage 
one A cell into a cancerous stage two cell; the division 
rate fi of the stage two cells; and the degree of imbalance 
A between their stochastic rate of proliferation and dif- 
ferentiation. In summary, focusing on the proliferating 
cell compartment only, and denoting the stage two mu- 
tated cells as type A*, then the revised cell proliferation 
model includes the additional non-equilibrium processes 

A ^ A* 

4 * m (A*+A* prob. (l + A)/2 (10) 
\ prob. (1 - A)/2. 

The rate v may be interpreted as the mean rate with 
which a stage-one cell acquires an additional mutation 
necessary to activate a second oncogene. The mutated 
cells then give rise, on average, to an exponentially grow- 
ing cell lineage with growth rate A/x. 

This nonequilibrium process was originally addressed 
by Kendall, who predicted the distribution in the num- 
ber of tumours detected at time t after mutation [35| . His 
focus on tumour statistics may reflect the experimental 
limitations in clonal analysis at the time: Until recently 
it was not possible to reliably detect clones at all, let 
alone to count the number of cells per clone. Experimen- 
tally, however, the clone size distributions are a more ef- 
ficient measure of cell kinetics than the tumour number 
distributions, because they result in a far richer data set, 



11 



and are accessible within weeks rather than months. We 
shall therefore extend Kendall's approach to predict the 
clone size distributions at times far earlier than tumour 
appearance. 



to the Master equation below, one may show that the 
variance in clone size also changes qualitatively: Whereas 
for normal skin the RMS variance in clone size grows as 
i 1 / 2 , here the variance in the long-time limit is finite, 



B. Clonal behaviour during early-stage cancer 

To familiarise ourselves with the modified model, con- 
sider the evolution of the average clone size with time. 
Focusing on the proliferating cell compartment with n 
type A cells and n* type A* cells in a clone, the relevant 
mean-field equations are 

d t {n) = -v{n) , 

dt(n*) = v(n) + A/x(n*) , 

which give the expected shift from linear growth of clones 
in normal skin to that of exponential growth, (n + n*} — 
(i/e A/it + A/ie _i/t )/(^ + A/i). More interestingly, referring 



lim 

t — >oo 



((n* - (n*)) 2 ) 1 / 2 



That is, the relative broadening of the clone size distribu- 
tion observed in normal skin is halted by the introduction 
of an exponentially growing cell population. 

These observations may already provide a crude 
method for identifying carcinogenesis through clonal 
analysis. To do better, it becomes necessary to solve for 
the full size distribution by extending the Master equa- 
tion (J3J) to include process (fTOjl . If we neglect the fate 
of differentiated cells, then the Master equation now de- 
scribes the evolution of the probability P n ,n* (t) for find- 
ing n type A cells and n* type A* cells in a clone, 



d t P n ,n* = r\ [(n - l)P n _i, n * - uP„,„«] + rX [(n + l)P n +i,n« - nPn.n*] + v [(n + l)P n+ i, n »_i - rcP„ >n « 
+ [( n * - l)Pi,n*-i - n*P n , n *] + /i [(n* + l)P„, n » + i - n*P n , n *} , 



r 



subject to the experimental boundary condition 
P«,n*(0) = <5i,o corresponding to exactly one "stage one" 
cell per clone at t = 0. As for the case of normal skin, we 
shall later be interested in the distribution of persistent 
clones, defined as, 
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P n ,N- n{t) 
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^ l-Poo(t) 



Pi,o(t) - P ,i(t) 



While it is not possible to solve Eq. (|11[) analytically, 
progress may be made when we allow for the widely- 
accepted view that tumours are monoclonal, that is they 
arise from a single "stage two" mutated cell 33]. This 
assumption conveniently limits us to the parameter space 
v <C A/it, for which an approximate long-time solution for 
the full clone size distribution may be found. 

Referring to the appendix for details, we find that the 
binned clone size distribution takes the long-time asymp- 
totic scaling form, 



P r=> w ^[i A .(^) 

where <j> k (t) = (1 + A- 1 )e At * t /2 k 
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Despite its ap- 



Kf - 4rA x 2 2 

13 = 2 X rX/A», and a = §^ 
parent complexity, this distribution is characterised 
by a simple scaling behaviour: Referring to Fig. EJa), 
the predicted clone size distributions are plotted using 
the scaling appropriate to the normal (unperturbed) 
system (cf. Fig. |4}. In this case, it is apparent that 
the scaling t i— > t/2 k fails. By contrast, from the 
expression for 4>k{t), it is clear that the size distribu- 
tions should scale according to the time translation, 
t i ► t' k = t + fcln2/A/x as confirmed by the results 
shown in Fig. [7(b). Further consideration of the size 
distribution exposes several additional features, which 
may provide further access to the new model parameters: 



• The long-time distribution decays with a rate /3A/x: 

Expanding Xg ja (x) for small x gives us the asymptotic 
form of the universal decay curve. For (3 < 1, consistent 
with the monoclonicity requirement A/i 3> v, we find 



lim V 

t>A/i 



(cane.) 



(12) 



where T(x) denotes the Gamma function. This expres- 
sion allows us to estimate (3 from the rescaled clone size 
distributions, providing access to the cell division and 
mutation parameters of the observed cells. 
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FIG. 7: (Color online) (a) The total number of basal layer cells per labelled clone during the onset of cancer according to 
process f|10p . The figure was plotted by numerically integrating Eq. (|ll|l using the empirical value rX = 0.088/week found 
for normal skin, together with hypothetical values of the cancer growth parameters v = 0.1 rX, /j, = 10 rX, and A = 0.5. To 
compare with normal skin, the predicted clone size distributions are replotted against the rescaled time coordinate t/2 k i— » t 
in (a) inset. In contrast with Fig. [4] here the curves fail to converge. In (b), the same curves are shown converge onto the 
universal form given in Eq. pip (dashed) when they are plotted against a new rescaled time 1 1— ► i! k — t + fcln2/A/^. Note that 
the large-clone distributions converge rapidly, whereas the distributions for smaller clones are affected by the non-negligible 
contribution of non-cancerous (A) cells to the small-clone size distribution. 



• The probability of tumour formation is finite: 

This is a well-known feature of the simple non-critical 
birth-death process (fTfJ|) [27} . Referring to the appendix, 
we find that the probability px for any given clone to 
survive and form a tumour is finite, 



PT = 1+ 27A" 



// v \2 2Ai/ 
\2r\) + r\(l + A) 



As a result, the onset of cancer will halt the steady de- 
crease in the density of labelled clones that is a hallmark 
of the unperturbed system. 

These properties, and especially the change in scal- 
ing behaviour, allow the onset of early-stage cancer to 
be identified from observations of clones less than one 
hundred cells in size. This may provide a dramatic im- 
provement both in speed and accuracy over current ex- 
perimental models, which rely on much later observations 
of tumours (or hyperplasias) in order to deduce the cell 
kinetics at early-stages. 



IV. CONCLUSIONS 

To summarize, we have shown that the range of clone 
fate data obtained from measurements of murine tail epi- 
dermis are consistent with a remarkably simple stochas- 
tic model of cell division and differentiation involving just 
one proliferating cell compartment. These findings over- 
turn a long-standing paradigm of epidermal fate which 



places emphasis on a stem cell supported epidermal pro- 
liferative unit. As well as providing significant new in- 
sight into the mechanism of epidermal homeostasis, these 
results suggest the utility of inducible genetic labelleling 
as a means to resolve the mechanism of cell fate in other 
tissue types, and as a means to explore quantitatively the 
effects of drug treatment and mutation. 

To conclude, we note that the analysis above has fo- 
cused on the dynamics of the clonal population with- 
out regard to the spatial characteristics. Indeed, we 
have implicitly assumed that any model capable of de- 
scribing the cell size distributions will also succeed in 
maintaining the near-uniform areal cell density observed 
in the basal layer. However, it is known that, when 
augmented by spatial diffusion, a simple Galton- Watson 
birth-death process leads to "cluster" formation in the 
two-dimensional system whereupon local cell densities di- 
verge logarithmically [36|, l37[ . Significantly, these diver- 
gences can not be regulated through a density-dependent 
mobility. Understanding how the Galton- Watson process 
emerges from a two-dimensional reaction-diffusion type 
process represents a significant future challenge. 

From a practical perspective, there is also the signif- 
icant question of how the cell kinetic model might be 
generalised to describe other forms of epidermis. In par- 
ticular, it is not feasible to repeat these experiments in 
vivo in human epidermis, a system of obvious medical 
significance. Therefore, it may be of great interest to de- 
termine, in future studies, the extent to which our results 
compare with the behaviour found in other systems. 

Lastly, our analysis of the cancer system referred to 
the relatively simple case of a two-stage mutation. It is, 



13 



of course, well-known that tumour formation is usually 
the result of multiple mutations. Understanding whether 
clonal fate data can be used to probe the kinetics of multi- 
stage mutation remains an interesting future challenge. 

Acknowledgements: We are grateful to Sam Ed- 
wards, Martin Evans, Bill Harris, Marc Kirschner, and 
Gunter Schiitz for useful discussions. The experimental 
aspect of this work was funded by the Medical Research 
Council, Association for International Cancer Research 
and Cancer Research UK. 



APPENDIX: CLONE SIZE DISTRIBUTIONS IN 
THE TWO-STAGE CANCER MODEL 



To derive the clone size distribution given in Eq. (|11| . 
we start by quoting the known result for the probability 
distribution n„» (t; r) of finding n* stage-two A* cells at 
time t starting from a single A* cell at time r [13] , 
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n„,(t;r) 



t) 3> 1, this distribution asymptotes to the 



( (!~ A ) 
IT+AJ 



( 



, 2A \ -An(t-T 
{JT+A)J e 

I * 2A 

exp(-n (T^sy. 



-Afi(t- 



for n* = 



for n* > 1 



(14) 

From the value of Llo we see that even when a cell has 
mutated, it is not guaranteed to result in a tumour: This 
will only occur with a probability of / = 1 — LTo(< — » 
oo) = (2A)/(1 + A). The value of / plays an important 
role in determining the statistics of tumour formation, as 
will be seen below. 

We now make two approximations: First, we take the 
long-time clone size distribution to be dominated by the 
statistics of A* cells. This is a safe assumption at times 
t > 1/u and t ^> 1/A/u, as may be seen by consid- 
ering the behaviour of the mean-field equations in sec- 
tion IIII Bl This approximation allows us to focus on 
the size distribution of A* cells only, p n - (t) , which is 
related to the full clone size distribution by the sum 
Pn*{t) = X)^Lo Pn,n* (*)• Secondly, we assume that the 
entire population of type A* cells in each clone arises 
from the first mutated cell that gives rise to a stable, ex- 
ponentially growing lineage of cells. This corresponds to 
the condition v -C A/i, as discussed in the main text. 

With these two approximations, the probability of find- 
ing a labelled clone containing n* > mutated cells is 
given by the population distribution of the first surviving 
cell lineage of A* cells, 

p n *(t)~N / TLn.it -t) V(l-/) m - V m (r)dr, (15) 
J ° m=l 

where J\f is some normalisation constant, and we have in- 
troduced the probability r m (r)dr for the m-th lineage of 



mutated cells within a given clone to be created during 
the interval r — > r + dr through the mutation process 
A — > A*. The weight factor (1 - f)" 1 ^ 1 gives the proba- 
bility that the first m — 1 cell lineages of A* cells within 
a clone will become extinct — a situation necessary to 
make the m-th cell line relevant to the distribution ac- 
cording to the monoclonal approximation. 

The rates r m (r) may be accessed by considering the 
probability w nim (t) that a clone containing n type A 
cells at time t also contains m independent lineages 
of mutated A* cells, each arising from a separate mu- 
tation event. (Later we shall treat the evolution of 
each of these cell lines post-creation). To solve for 
Wn m(t) we must introduce its moment- generating func- 
tion G(q,Q*;t) = £„=o Em=o «>n,m(t)« n (Q*) m , which 
evolves (from Eq. Hip according to the dynamical equa- 
tion, 



G= [r\{q-lf + u{Q* -q)]d q G. 



(16) 



Solving this equation subject to the initial condition of 
one "stage-one" (A) cell per clone, we find the solution 



G(q,Q*;t) = fo 



«Q*-(9- 1 -27a) 
+ 1 + 7^7, (17) 



with f g . ee y J{^ + JL.{i-Q*). Eq. (T7|) describes 
the evolution of a single A cell as it proliferates and even- 
tually gives rise to a set of internal lines of mutated cells. 

Before we proceed to find p n * (t), note that setting q — 
1, Q* = 1 — / in Eq. (fT7|) gives us the result (quoted in 
the main text) for the asymptotic fraction pt of clones 
in which all mutated cell lines become extinct, 



PT = 1+ 2r~\- 



\2rXJ rX 
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On the other hand, setting q — 1 only in Eq. (fTT)) 
gives the moment-generating function for (yet another) 
distribution W m (t) = ^2 n=0 w ntm (t) of a clone containing 
?7i independent lines of A* cells irrespective of the number 
of normal cells in the clone. Finally, noting that W m (t) = 
r m - r m+ i then gives: 

m— 1 

r m (t) = - W m (t), 



which we may substitute into Eq. (|15[) to find (for n* > 0) 

p n *(t)~~JV f drn„.(t-r)G(l,l-/;r). (18) 
Jo 

From this expression, simplified by the large-n* approx- 
imation (Y) n » — J dn), we obtain the final form of the 
binned size distribution given in Eq. (jll[) . 
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